##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--muti_cancer"), type = "character") ,
    make_option(c("--muti_pre"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    muti_cancer <- "~/20220915_gastric_multiple/dna_combinePublic/maf/All_GGA.cancer.maf"
    muti_pre <- "~/20220915_gastric_multiple/dna_combinePublic/maf/All_GGA.precancer.maf"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/ITH"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
muti_cancer <- opt$muti_cancer
muti_pre <- opt$muti_pre
sample_info <- opt$sample_info
out_path <- opt$out_path

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_mutiCancer <- fread(muti_cancer  ,sep = "\t",  quote="" ,header = T)
dat_mutiPre <- fread(muti_pre  ,sep = "\t",  quote="" ,header = T)
dat_mutiNodule <- rbind(dat_mutiCancer , dat_mutiPre) 

dat_info <- data.frame(fread(sample_info))

###########################################################################################

col<-brewer.pal(12,"Set3")
col_point <- brewer.pal(9,"Set1")
Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
dat_mutiNodule <- subset(dat_mutiNodule , t_alt_count > 0)

dat_mutiNodule <- data.frame(Hugo_Symbol = dat_mutiNodule$Hugo_Symbol,
    Chromosome = dat_mutiNodule$Chromosome , Start_Position =  dat_mutiNodule$Start_Position , End_Position = dat_mutiNodule$End_Position ,
    Reference_Allele = dat_mutiNodule$Reference_Allele , Tumor_Seq_Allele2 = dat_mutiNodule$Tumor_Seq_Allele2 , 
    Variant_Classification = dat_mutiNodule$Variant_Classification,
    t_ref_count = dat_mutiNodule$t_ref_count , t_alt_count = dat_mutiNodule$t_alt_count ,
    Tumor_Sample_Barcode = dat_mutiNodule$Tumor_Sample_Barcode)

dat_mutiNodule$Location <- paste( dat_mutiNodule$Chromosome , dat_mutiNodule$Start_Position , 
    dat_mutiNodule$Reference_Allele , dat_mutiNodule$Tumor_Seq_Allele2 , sep=":"  )

###########################################################################################
## 计算每个Normal的Tumor的瘤内异质性
## 共享的突变数目/各自私有的突变数目之和
## Genomic comparison of esophageal squamous cell carcinoma and its precursor lesions by multi-region whole-exome sequencing
## 2017. Nature Communications
###########################################
## Function
computeITH <- function( tmp1 , tmp2 ) {
    shareMut_Num <- length(which(tmp1$Location %in% tmp2$Location))
    t1_private_Num <- nrow(tmp1[!(tmp1$Location %in% tmp2$Location),])
    t2_private_Num <- nrow(tmp2[!(tmp2$Location %in% tmp1$Location),])

    ITH <- 1 - shareMut_Num/(t1_private_Num + t2_private_Num + shareMut_Num)

    res_tmp <- data.frame( shareMut_Num = shareMut_Num , t1_private_Num = t1_private_Num , t2_private_Num = t2_private_Num , ITH = ITH )

    return(res_tmp)
}

###########################################
resulst_hetetogeneity <- c()

for( Sample in unique(dat_info$ID) ){
    tmp <- subset( dat_info , ID == Sample )

    ## 计算ITH
    for( i in 1:(dim(tmp)[1]-1) ){
        for( j in (i+1):dim(tmp)[1] ){
            tumor1 <- tmp[i , "Tumor"]
            tumor2 <- tmp[j , "Tumor"]

            class1 <- tmp[i , "Class"]
            class2 <- tmp[j , "Class"]

            tmp1 <- subset( dat_mutiNodule , Tumor_Sample_Barcode == tumor1 )
            tmp2 <- subset( dat_mutiNodule , Tumor_Sample_Barcode == tumor2 )

            res_tmp <- computeITH( tmp1 , tmp2 )

            res_tmp2 <- cbind(data.frame( Normal = Sample , Tumor_1 = tumor1 , Tumor_2 = tumor2 ,Class_1 = class1 , Class_2 = class2) , res_tmp )

            resulst_hetetogeneity <- rbind( resulst_hetetogeneity , res_tmp2 )
        }
    }
}

out_name <- paste0(out_path , "/ITH.compute.allSample.tsv")
write.table( resulst_hetetogeneity , out_name , row.names = F , sep = "\t" , quote = F )

###########################################
## https://molecular-cancer.biomedcentral.com/articles/10.1186/s12943-019-0939-9#MOESM5
## 每个患者得到唯一的瘤内异质性索引，所有的ITH加起来/组合的数量
## IM + IGC + DGC的样本排除分析
## 只考虑IM_GC间的
exclusion_sample <- unique(subset( dat_info , Type=="IM + IGC + DGC")$ID)
result_hetetogeneity_unique <- subset( resulst_hetetogeneity , Class_2 != Class_1 & !( Normal %in% exclusion_sample ) )

result <- c()
for( id in unique(result_hetetogeneity_unique$Normal) ){

    tmp <- subset(result_hetetogeneity_unique , Normal == id )
    sample_num <- length(unique( c(tmp$Tumor_1 , tmp$Tumor_2 )))

    n <- sample_num
    r <- 2
    ## 用于计算从sample_num个元素中选择2个元素的组合数：
    combination <- factorial(n) / (factorial(r) * factorial(n - r))
    combination <- nrow(tmp)
    ith_sum <- sum(tmp$ITH)

    tmp <- data.frame( ID = id , ITH = ith_sum/combination )
    result <- rbind( result , tmp )
}

result_hetetogeneity_unique <- result
result_hetetogeneity_unique <- merge( result_hetetogeneity_unique , unique(dat_info[,c("ID" , "Type")]) )

## IM + IGC + DGC的样本算两份
result_hetetogeneity_igc_dgc <- subset( resulst_hetetogeneity , Class_2 != Class_1 & ( Normal %in% exclusion_sample ) )
result <- c()
for( id in unique(result_hetetogeneity_igc_dgc$Normal) ){

    ##################
    tmp <- subset(result_hetetogeneity_igc_dgc , Normal == id & Class_1 != "DGC" & Class_2 != "DGC" )
    sample_num <- length(unique( c(tmp$Tumor_1 , tmp$Tumor_2 )))

    n <- sample_num
    r <- 2
    ## 用于计算从sample_num个元素中选择2个元素的组合数：
    combination <- factorial(n) / (factorial(r) * factorial(n - r))
    combination <- nrow(tmp)
    ith_sum <- sum(tmp$ITH)

    tmp <- data.frame( ID = id , ITH = ith_sum/combination )
    tmp$Type <- "IM + IGC + DGC(IGC)"
    result <- rbind( result , tmp )

    ##################
    tmp <- subset(result_hetetogeneity_igc_dgc , Normal == id & Class_1 != "IGC" & Class_2 != "IGC" )
    sample_num <- length(unique( c(tmp$Tumor_1 , tmp$Tumor_2 )))

    n <- sample_num
    r <- 2
    ## 用于计算从sample_num个元素中选择2个元素的组合数：
    combination <- factorial(n) / (factorial(r) * factorial(n - r))
    combination <- nrow(tmp)
    ith_sum <- sum(tmp$ITH)

    tmp <- data.frame( ID = id , ITH = ith_sum/combination )
    tmp$Type <- "IM + IGC + DGC(DGC)"
    result <- rbind( result , tmp )
}

result_hetetogeneity_unique_final <- rbind( result_hetetogeneity_unique , result )
result_hetetogeneity_unique <- result_hetetogeneity_unique_final

out_name <- paste0(out_path , "/ITH.compute.uniqueNormal.tsv")
write.table( result_hetetogeneity_unique , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 比较IM和IGC还有DGC的瘤内异质性的差异

p <- wilcox.test(subset(result_hetetogeneity_unique , Type=="IM + IGC")$ITH , subset(result_hetetogeneity_unique , Type=="IM + DGC")$ITH)$p.value

if( p < 0.001 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
}

result_hetetogeneity_unique$p_text <- p_text

plot <- ggplot(data=result_hetetogeneity_unique,mapping = aes(x=Type,y=ITH))+
  geom_boxplot(lwd=1.5,aes(color=Type) , outlier.colour = NA) +
  geom_jitter(position=position_jitter(0.2),aes(color=Type)) +
  scale_color_npg() +
  #facet_grid(.~Type)+
  xlab(NULL) +
  ylab('ITH')+
  theme_bw() +
  ylim(0.2,1) +
  geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
  #stat_compare_means(label.y = 1.1) +
  theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='none',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        legend.text = element_text(size = 8,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 12,color="black",face='bold'),
        axis.line = element_line(size = 0.5)) 

out_name <- paste0(out_path , "/ITH.combineNormal.IGC_DGC.pdf")
ggsave(file=out_name,plot=plot,width=3,height=5)